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Abstract 


A  simple  form  of  consistency  relations  between  generalized  forces  and  displacements 
for  systems  exhibiting  power-law  behavior  is  presented.  The  later  discussion  focuses 
on  certain  details  regarding  applications  of  the  relations  to  nonlinear  fracture  mechan¬ 
ics,  emphasizing  the  finite  element  analysis  of  a  single  edge-cracked  strip  subjected  to 
remote  tension  under  plane  strain  conditions.  - - 
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Introduction 


Within  the  last  decade  there  has  been  observed  a  serious  computational  effort  to¬ 
wards  obtaining  numerical  solutions  to  problems  of  nonlinear  fracture  mechanics.  The 
handbook  formatted  reference  (Kumar  et  al.,  1981^  gives  extensive  tabulations  of  im¬ 
portant  parameters  (J-integral,  CTOD,  eb^for  various  configurations,  loadings  and 
material  behaviors,  which  makes  it  useful  for  engineering  analysis,  and  the  review  of 
prior  work  in  the  field  permits  its  use  it  as  a  good  starting  point  for  further  research. 
The  later  work  of  Shih  and  Needleman  (1984)  exposes  some  quantitative  discrepancies 
between  the  two  results.  These  discrepancies  cause  confusion  and  implicitly  raise  the 
legitimate  question:  ‘Who  is  right?*.  The  purpose  of  this  work  is  to  help  to  resolve  the 
matter  of  quality  of  reported  data. 

There  are  three  parts  in  this  paper.  The  first  one  focuses  on  the  derivation  of 
consistency  relations  between  generalized  forces  and  displacements  in  fairly  general 
types  of  nonlinear  systems,  including  the  class  of  traction  prescribed  boundary- value 
problems  for  power-law  isotropic  materials,  exhibiting  in  pure  tension  stress-strain 
behavior 


where  to  and  <To  ve  reference  strain  and  stress  values,  and  n  is  a  material  exponent 
varying  from  one  for  the  linear  material  to  infinity  for  the  rigid  plastic  material. 

The  second  section  discusses  some  specific  features  of  solutions  for  fracture  me¬ 
chanics  configurations,  stressing  both  the  validity  of  the  application  of  the  consistency 
relations  in  principle  and  their  potential  use  in  the  detection  of  errors. 

The  paper  is  concluded  with  a  numerical  example  to  which  the  consistency  con¬ 
ditions  are  applied  -  the  single  edge<racked  strip  subjected  to  remote  tension  under 
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plane  strain  conditions.  The  results  are  compared  against  the  data  reported  by  Kumar 
et  al.  (1981)  and  Shih  euid  Needleman  (1984). 

Throughout  the  paper  we  use  boldface  letters  for  vectors,  matrices  and  tensors. 
Dot  stands  for  an  appropriate  inner  product.  All  computations  reported  herein  have 
been  performed  using  the  ABAQUS  finite  element  program  on  a  DATA  GENERAL 
MV- 10000  computer. 


Consistency  Relations 


Let  us  consider  some  body  loaded  by  a  set  of  generalized  forces  Q.  The  set  of 
generalized  displacements  is  taken  as  derivable  from  the  constitutive  potential  F,  which 
depends  on  the  overall  geometry,  material  properties  and  the  forces: 


q  = 


dQ' 


(2) 


We  confine  our  attention  to  the  class  of  potentials  which  are  both  convex  and  homo¬ 
geneous  functions  of  Q  degree  (n  4- 1): 


F  (aQi  +  (1  -  a)  Q,)  <  aF  (Qi)  +  (1  -  a)  (Qa) , 
||.Q  =  (n^l)F. 


(3) 

(4) 


In  equation  (3)  it  is  understood  that  Qi.Qa  are  arbitrary  sets  of  generalized  forces 
and  0  <  o  <  1.  In  the  relation  (4),  we  adopt  Euler’s  Theorem  on  homogeneous 
functions  aa  a  definition.  One  of  the  important  properties  of  convex  functions  is  the 
positive  semidefiniteness  of  the  Hessian  matrix  T  (tangent  compliance  in  the  force- 
displacement  framework): 


r  = 


d^F 

dQ^' 


(5) 
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The  inner  product  of  F  with  Q  ®  Q  leads  to  an  important  conclusion  that  the  potential 
functions  F  defined  above  are  nonnegative. 

The  adternative  definition  of  a  homogeneous  function  allows  us  to  reduce  the  effective 
number  of  independent  variables  by  one;  the  most  general  representation  serving  our 
purposes  is: 

f  (Q)  =  /(x)  (r(q)l”* ,  (6) 


where  /(x)  is  a  nonnegative  function  of  the  reduced  set  x,  which  dimension  is  one 
less  than  that  of  Q,  and  r(Q)  is  an  arbitrarily  defined  norm  in  the  Q-space.  The  key 
point  of  the  forthcoming  derivation  is  that  the  generalized  displacements  in  equation 
(2)  are  the  partial  derivatives  of  the  same  function  amd,  therefore,  should  be  related. 
We  limit  our  demonstration  to  the  case  of  two  generalized  forces.  The  extension  to 
a  multidimensional  set  is  transparent  and  will  not  be  considered  here.  In  this  case 
equation  (6)  may  be  reduced  to 

9:=<3.>0,  x:=|i.  (7) 

Straightforward  differentiation  leads  to  the  following  expressions  for  q  and  F: 


I  r  )n^V 


(8) 


/  (n  +  l)n/  -  2n*/'  +  x’/"  symmetry 

\  nr-xr  r  jn  +  r 


(9) 


where  prime  denotes  differentiation  with  respect  to  x.  To  enforce  positive  semidef¬ 
initeness  of  F  we  require  both  Fj]  and  determinant  of  F  to  be  nonnegative,  which 
gives 


r  >0, 


(10) 
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ir  -  >  0. 

n  -t- 1 


(11) 


But  since  the  function  /(x)  is  nonnegative,  condition  (11)  is  sufficient  to  guarantee 
positive  semidefiniteness. 

Introduction  of  an  auxilary  function  y,  in  a  way  a  conjugate  of  x,  and  given  by 


9a 

y  :=  —  = 


r 


qi  (n  +  l)/-x/'’ 
permits  the  rewriting  of  condition  (11)  in  a  simpler  form: 

y'>0. 


(12) 


(13) 


Now,  instead  of  considering  the  related  pair  of  functions  /(x)  and  /'(x)  to  describe 
the  generalized  displacements,  we  can  use  the  pair  /(x)  and  y(x)  related  via  (12).  If 
values  of  /(x)  and  y(x)  are  given  at  some  point,  say  x  =  Xo,/(xo)  =  /o,y(®o)  =  yo> 
then  integration  of  (12)  leads  to 

/(x)  =  /„«,((n  +  l)/Jj^).  (H) 

But  for  any  value  of  x  >  xo  with  xq  <  t  <  x  we  may  write  the  following  inequality 


yo  <  y(0  <  y{^) 

1  +  tyo  ~  1  +  ty{t)  “  1  +  ty{x)  ’ 


(15) 


which,  after  the  integration,  gives  us  the  main  result  of  this  section: 

(16) 

1  +  xolto  \fo)  1  +  xoy(x) 

It  is  important  to  mention  that  (11)  and  (13)  are  equivalent  only  if  yi  >  0  for  all 
the  points  within  the  interval  [xo,x],  which  implicitly  imposes  conditions  on  both  xo 
and  X. 
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Application  to  Fracture  Mechanics 


A  broad  class  of  problems  of  fracture  mechanics  can  be  described  by  a  generic 
problem  -  one  is  given  a  configuration  containing  a  crack,  loaded  by  forces  Q,  and  is 
asked  to  determine  the  value  of  the  J-integral  (Rice,  1968)  to  characterize  the  local 
fields  (Hutchinson,  1968,  Rice  and  Rosengren,  1968)  and  q  -  displacements,  at  the 
remote  distances,  due  to  the  introduction  of  the  crack, 

q  =  q‘  -  q\  (17) 

where  q‘  and  stand  for  generalized  displacements  in  two  auxilary  problems.  The 
first  one  is  a  prescribed  traction  boundary-value  problem  for  the  given  configuration. 
The  second  one  is  identical,  but  there  is  no  crack. 

The  constitutive  potential  F(Q)  becomes  equal  to  the  difference  in  the  complemen¬ 
tary  energies  of  the  auxilary  problems: 

F(Q)  =  -L-  f  (E‘ .  T*  -  E*  •  T*)dV,  (18) 

n  +  1 

where  T  and  E  stand  for  the  stress  and  strain  tensors,  respectively. 

The  connection  between  this  class  of  problems  and  the  one  described  in  the  first 
section  is  obvious  but,  nevertheless,  there  are  some  important  details  to  be  considered 
both  due  to  the  necessity  to  perform  numerical  analysis  and  the  special  features  of  the 
fractxire  mechanics  problems  per  $e. 

At  first,  we  would  like  to  address  the  matter  of  possible  loss  of  homogeneous  struc¬ 
ture  of  F{Q).  In  principle,  the  linearity  of  equilibrium,  compatibility  equations,  and 
boundary  conditions,  combined  with  the  constitutive  law,  requires  an  analytical  solu¬ 
tion  to  the  problem  to  be  homogeneous.  But  if  we  have  to  model  an  incompressible 
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material  and,  therefore,  employ  a  penalty  procedure  in  numerical  analysis,  we  can  en¬ 
counter  deviations  from  homogeneity.  The  simplest  way  to  introduce  a  penalty  is  by 
linear  relations  between  isotropic  components  of  the  stress  and  strain  tensors: 

trT  =  ^trE.  (19) 

3co 

We  argue,  hueristically  at  best,  that  as  K  tends  to  infinity,  the  infiuence  of  the  hy¬ 
drostatic  stress  diminishes,  and  material  response  tends  toward  incompressibility.  An 
attempt  to  conserve  homogeneity  by  using  a  power-law  penalty  may  easily  lead  to 
numerical  problems  as  we  operate  with  large  numbers.  Therefore,  we  can  claim  that 
degree  of  compressibility  and  deviation  from  the  homogeneous  structure  are  implicitly 
related,  and  in  the  limiting  case  of  a  large  material  tends  towards  both  homoge¬ 
neous  and  incompressible  response.  The  straightforward  application  of  (16)  may  be  of 
use  in  detecting  errors  in  interpretation  of  penalty  term  for  sophisticated  variational 
formulations  combining  both  regular  displacement  based  and  hybrid  (displacement  and 
pressure)  formulations  of  the  finite  element  method,  as  implemented,  for  example,  in 
ABAQU S.  Later  on  we  present  a  numerical  example  of  this  relation. 

The  question  of  convexity  is  especially  interesting  for  the  fracture  mechanics  prob¬ 
lems.  The  difficulty  here  is  that  there  is  no  single  boimdary-value  problem  which  may 
be  directly  analyzed  to  determine  /'(Q),  but  rather  two  auxilary  problems.  Convexity 
of  each  problem  is  guaranteed  by  the  monotonicity  of  the  stress-strain  curve  or,  more 
precisely,  by  the  convexity  of  the  strain  energy  density  function  (Marsden  and  Hughes, 
1983).  The  function  F  is  nonnegative  (Rice,  1968),  but  is  not  necessarily  convex.  The 
straightforward  mathematical  example  is  F(Qi,Qs)  =  |QiQ3|°^- 

In  order  to  demonstrate  the  possible  loss  of  convexity  in  a  physical  problem,  we 
address  the  case  of  a  penny-shaped  crack  embedded  in  an  infinite  isotropic  power- law 


6 


matrix  subjected  to  axisymmetric  remote  loading,  characterized  by  the  axial  and  radial 
stresses  Qi  and  respectively.  We  consider  three  materials  which  can  be  generated 
from  (l)  as  tensorial  extensions.  The  first  one  is  the  incompressible  material  given  by 


where  S  is  the  stress  deviator  and  we  take  n  =  3.  The  second  material  is  linear  elastic 
(n=l}  with  Yoimg’s  modulus  E  =  rro/eo  Poisson’s  ratio  u.  The  last  case  is  an 
isotropic,  power-law  material,  comprtsaibU  material,  which  constitutive  equations  are 
derived  from  (20), (21)  by  substitution  of  the  tensor  T  itself  for  the  stress  deviator 
S.  The  values  of  <ro,eo  ^d  u  are  not  important  as  they  are  appearing  as  constant 
multipliers. 

The  first  example  is  treated  in  detail  by  He  and  Hutchinson  (1981),  (1983);  for  the 
case  Q\  >  0,Qi  >  Q%,  The  suggested  functional  fit  to  the  numerical  solution  is  given 
by 

where  a  is  the  crack  radius.  The  above  repression  gives  us  concave  F(Q),  which  can 
be  seen  directly  from  deriving  the  Hessian  matrix.  The  formula  is  derived  from  a 
perturbation  technique,  agrees  well  with  numerical  solutions  providing  x  <  0.6,  and 
for  larger  x  it  fails  to  give  an  accurate  estimate  (He  and  Hutchinson,  1983).  The  finite 
element  analysis  of  the  problem,  which  we  have  conducted,  shows  that  at  approximately 
the  same  point  (z  »  .6)  the  complete  numerical  solution  gains  convexity;  therefore,  we 
may  conclude  that  there  is  conditional  convexity  in  this  case. 
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The  well-known  solution  (Sneddon,  1964)  for  the  isotropic  linear  material  coincides 
with  (22)  for  n  =1  and  u  -  -  and  is  given  by 

m)  =  ‘■•q!-  (23) 

This  expression  gives  only  one  non-zero  component,  Tn,  of  the  Hessian  matrix,  which, 
of  course,  retains  positive  semidefiniteness,  but  relations  (16)  degenerate  to  triviality. 

By  performing  a  finite  element  analysis  we  find  that  the  third  material  gives  us  an 
unconditionally  convex  potential,  and  the  main  result  of  the  first  section  is  of  use. 

In  the  above  examples  we  have  encountered  thr' s  possible  situations  for  two-dimensional 
Q-space;  namely,  relations  (16)  are  relevant  at  some  points  in  the  domain,  they  are 
relevant  throughout  Q-space,  or  they  can  not  be  applied  at  all.  The  heuristic  conclu¬ 
sion  may  be  drawn  if  we  consider  the  dimensions  of  the  Q-space  of  the  two  auxilary 
and  the  main  problem  for  all  three  materials.  In  the  first  example  the  solution  of  the 
problem  without  the  crack  depends  only  on  the  applied  equivalent  Mises  stress  &.  If 
we  introduce  a  new  set  of  generalized  forces,  namely 

Qi  =  =  9l  - (24) 

where  the  first  equation  of  (24)  defines  applied  hydrostatic  pressure  and  the  second 
the  equivalent  Mises  stress, then  the  dimension  of  Q'-space  is  one  in  the  context  of 
the  problem  without  the  crack.  The  main  and  the  other  auxilary  problem,  on  the 
other  hand,  exhibit  two-dimensional  load  space.  It  is  clear  &om  (23)  that  an  analogous 
situation  occurs  with  the  linear  elastic  material,  except  that  the  main  problem  is  the 
one  which  has  the  reduced  space.  The  potential  due  to  the  introduction  of  »he  crack 
does  depend  on  the  single  force  Qi,  but,  by  superposition,  is  independent  of  Q3.  Only 
for  the  last  example  is  there  a  truly  two-dimensional  Q-space  for  all  three  problems. 
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It  is  clear  that  in  the  first  two  examples,  there  exists  some  set  Q,  which  essentially 
forms  the  null'Space  of  ^(Q)  for  the  main  or  the  auxilary  problem.  Obviously,  it  is 
hydrostatic  pressure  for  the  incompressible  material  and  stressing  parallel  to  the  crack 
plane  for  the  main  problem  in  the  linear  case.  The  loosely  defined  induction  is,  that 
in  order  to  employ  (16),  all  three  problems  must  have  the  two-dimensional  Q-space.  A 
more  rigorous  statement  would,  perhaps,  require  the  definition  of  F  on  the  complement 
of  the  null-space. 

The  conclm  'on  is  hueristic  but,  nevertheless,  seems  to  rehabilitate  the  'misbehavior’ 
of  the  otherwise  mathematically  ‘loyal’  equations. 

It  Is  important  to  note  that  this  conclusion  does  not  put  any  question  marks  on 
the  substantial  body  of  theoretical  and  computational  effort  (Budiansky  et  al.  1981, 
Rafalsky,  1985)  in  terms  of  the  derivations  and  the  implementation  of  a  variational 
principle  for  determination  of  F  directly,  because  here  we  deal  only  with  the  generalized 
force  space  of  boimdary  tractions,  rather  than  with  a  function  space  of  Ritz  procedure. 

To  conclude  this  section  we  would  like  to  mention  the  possibility  of  including  the 
J-integral  into  a  gradient  structure  analogous  to  (2),  which  leads  to  the  correlation 
between  near  and  far  fields  and,  of  course,  to  another  group  of  consistency  relations. 
The  procedure  that  was  initially  suggested  by  Parks  et  al.  (1983)  and  later  was  applied 
by  Shih  and  Needleman  (1984),  is  essentially  based  on  finite  difference  approximation 
of  the  gradient  scheme  for  planar  and  axysimmetric  problems. 

Numerical  Example 

A  single  edge-cracked  strip.  Figure  1,  subjected  to  a  remote  loading  under  plane 
strain  conditions  is  considered.  The  dimensions  are: 

a  10,  6  =  20,  L  =  60. 
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The  constitutive  behavior  is  modeled  by  (19)-(21)  with  numerical  values  of  material 
constants: 

<7o  =  1,  Co  =  li  K  —  10*,  ft  =  5. 

The  values  of  <7o  and  Cq  do  not  correspond  to  any  real  material,  but  as  long  as  we 
cire  able  to  scale  solutions,  the  idea  of  operating  with  more  computationally  convenient 
numbers  is  rather  helpful.  The  finite  element  discretization  is  given  in  Figure  2.  The 
plane  strain  eight-node  hybrid  element  with  nine  integration  points  and  a  bi-linear 
interpolation  for  pressure  is  employed. 

The  remote  forces  per  tinit  thickness  are  tensile  load  N  =  20  and  varying  positive 
bending  moment  M,  tending  to  close  the  crack.  We  are  interested  in  the  pure  tension 
solution  and  apply  bending  only  for  the  purx>oee  of  simulation  of  conditions  (16).  The 
generalized  forces  are  identified  as  Qi  :=  Af,  Qi  :=  M/6  with  corresponding  pair  of 
generalized  displacements  qi  :=  S  and  qt  :=  0b.  It  is  understood  that  the  generalized 
displacements  refer  to  the  contribution  due  to  the  crack,  as  discussed  in  the  previous 
section. 

The  potential  F{Q)  and  the  7-integral  are  taken  in  convenient  dimensional  forms 
F(Q)  =  —cM)  (-)  ,  (25) 

J  =  ao£obh{x)  •  (26) 

The  results  of  the  computations  are  given  in  ro6/e  1.  We  conclude  that  the  current 
analysis  gives  results  close  to  those  of  Shih  and  Needleman  (1984)  for  the  values  of 
the  rotation  and  the  7-integral;  the  displacement  is  somewhere  between  the  earlier 
reported  data.  It  is  worth  mentioning  that  the  main  point  of  discrepancy  between  the 
previous  data  is  in  the  near  field  quantities;  therefore,  our  results  rather  support  the 
data  reported  by  Shih  and  Needleman  (1984). 
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The  next  part  of  our  analysis  is  concerned  with  the  simulation  of  data  for  relations 
(16).  It  is  obvious  that  in  this  example  problem,  consistency  relations  can  be  of  impor¬ 
tance.  But  on  the  other  hand,  if  values  of  x  and  xo  are  sufficiently  close  we  expect  the 
interference  of  the  numerical  noise  in  the  tabulated  data  to  be  substantial;  therefore, 
we  have  to  decide  on  a  minimum  value,  z^,  such  that  for  all  x  —  xq>  Zm>i»  consistency 
relations  give  us  reasonable  conclusions. 

We  adopt  a  very  simplified  estimate.  The  difference  between  bounding  terms  in 
(15)  at  the  pure  tension  limit  zo  =  0  is  (z  —  Zo)(y  -  ys),  and  this  should  remain  positive 
in  the  moet  unfavorable  case.  As  we  have  four  kinematical  data  entries  in  the  above 
formula,  then  we  require 

y  -  Sfo  >  4eq,  (27) 

where  Sq  is  the  relative  error  in  q.  From  the  homogeneity  and  dimensional  considera¬ 
tions  we  can  take: 

«q  =  (28) 

where  sr  is  the  relative  error  in  nodal  reaction  force  in  the  virtual  work  sense,  and  in 
the  analysis  with  the  ABAQUS  program,  the  maximum  value  of  this  error  has  been 
set  at  10~*  .  In  expression  (28),  we  have  neglected  a  dimensionless  constant  multiplier 
expected  to  be  of  order  unity.  Relations  (27), (28)  form  the  implicit  conditions  on  Zmin- 
The  results  of  computations  are  given  in  Table  2  and  give  z^tn  =  8.10~^,  which  is 
in  a  fair  agreement  with  conditions  (27), (28)  which  give  the  value  of  z^un  =  3.10'^. 
Therefore,  we  may  claim  that  our  data  is  proIxUtly  acceptable. 

The  section  is  concluded  by  the  consideration  of  the  penalty  term  as  a  possibility 
for  the  loss  of  homogeneity.  In  this  set  of  computations  we  take  K  varying  from  1  to 
10*,  keeping  the  rest  of  the  material  constants  to  be  the  same. 
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The  results  of  Table  3  show  that  the  maximum  difference  for  various  K  >  100,  in 
Xmin,  •/-integral  zmd  q,  only  appears  in  the  fourth  digit.  These  results  suggest  that  one 
might  be  able  to  use  moderate  values  of  the  penalty  K.  There  are  clear  advantages 
to  such  a  procedure,  because  the  computations  converge  more  rapidly  and  there  is  less 
chance  to  encounter  numerical  difficulties.  For  small  values  of  K  (below  100)  the  most 
sensitive  parameter  turns  out  to  be  Xmim  though  the  physical  quantities  remain  within 
an  acceptable  variation. 
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Reference 

/ 

y 

h 

83.90 

-2.393 

309.9 

113.1 

-2.430 

407.4 

current 

99.00 

-2.677 

390.2 

Table  1.  Pure  Tension  Data,  n  =  5,K  =  10®,o/6  =  1/2. 


X 

/ 

y 

Lower 

Bound 

Intermediate 

Value 

Upper 

Bound 

.0000 

99.039237 

-2.677448 

- 

- 

- 

98.881223 

-2.676500 

1.000268 

1.000266 

1.000268 

98.720673 

-2.675612 

1.000535 

1.000537 

1.000536 

1^^ 

98.564139 

-2.674622 

1.000803 

1.000802 

1.000804 

.0004 

98.405337 

-2.673701 

1.001071 

1.001071 

1.001072 

.0005 

98.247815 

-2.672751 

1.001338 

1.001338 

1.001341 

.0006 

98.090018 

-2.671816 

1.001606 

1.001606 

1.001609 

.0007 

97.929625 

-2.670958 

1.001873 

1.001880 

1.001878 

.0008 

97.776636 

-2.669905 

1.002140 

1.002141 

1.002147 

97.618779 

-2.668991 

1.002408 

1.002411 

1.002416 

.0010 

97.461915 

-2.668057 

1.002675 

^002679 

1.002685 

Table  2.  Combined  Loading  Data,  K  =  10®,  n  =  5,o/6  =  1/2.^ 


logK 

h 

/ 

y 

®m*f» 

0 

393.7 

101.8 

-2.644 

1 

390.5 

99.30 

-2.674 

.0012 

2 

390.3 

99.05 

-2.676 

.0008 

4 

390.2 

99.05 

-2.677 

.0008 

6 

390.2 

99.04 

-2.677 

Table  3.  Pure  Tension  Data  for  Variable  Penalty  Term,  x  =  OjTI  =  b,afb  =  1/2. 

^  In  the  upper  block  of  the  table,  relations  (16)  are  not  valid,  in  the  next  block 
they  should  be  valid  according  to  (27), (28),  and  in  the  last  part  relations  (16)  are 
unconditionally  valid. 


